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A  Method  for  Computing  Spectral  Reflectance 


1.  Introduction.  The  colour  of  an  object  depends  on  the  spectral  distributions  of  the 
incident  light  and  the  reflectance  functions  of  the  object.  More  precisely  at  each  point  x 
the  receptors  (for  example,  the  rods  and  cone  of  the  eye)  measure  a  quantity  giver  by 

=  J  a^(\)i:(x,\)S(x,\)dk  (U) 

where  a„(X)  is  the  absorption  function  of  the  receptors  (n  =  I  to  4  if  we  consider  the  rods 
and  cones),  E(x,\)  is  the  incident  illumination  and  S(x,  X)  the  reflectance  of  the  surface. 
In  terms  of  the  notation  used  by  Hurlbert  (1905)  our  function  S(x,  X)  is  the  product  of  the 
albedo  p(x,  X)  and  the  reflectance  function  R(n,k,s). 

These  quantities  !„(*)  are  then  evaluated  and  combined  to  give  the  perception  of  colour. 
For  the  human  eye  it  is  assumed  that  the  outputs  of  the  three  cones  correspond  to  the 
perception  of  colour.  The  three  cones  corresponds  to  the  three  dimensions  of  the  space  of 
perceived  colours  (von  Helmholtz  1866,  or  see  Feynman  1965  for  an  interesting  review).1 

A  number  of  experiments  show  that  perceived  colour  of  an  object  is  relatively  independent 
of  the  spectral  distribution  of  the  incident  light  .  This  is  the  colour  constancy  effect 
demonstrated  most  clearly  by  Land  on  Mondrians.  Since  the  incident  illumination  and  the 
surface  reflectance  are  confounded  in  the  input  (1.1)  some  assumptions  must  be  made  to 
disentangle  them  and  to  obtain  a  colour  which  depends  on  the  surface  reflectance  alone. 
This  paper  suggests  a  scheme  for  doing  this.  Alternative  schemes  are  described  in  Rubin 
and  Richards  (1984),  Hurlbert  and  Poggio  (1985). 

This  paper  assumes  that  both  the  colour,  the  incident  illumination  and  the  surface  reflectance 
can  be  expressed  as  a  finite  sum  of  basis  functions  in  spectral  space.2  It  is  then  shown 
that  given  a  sufficient  number  of  patches  of  different  surface  reflectivity  with  roughly  the 
same  incident  illumination  it  is  possible  to  solve  for  the  surface  reflectance  up  to  an  overall 
scaling  factor.  In  such  a  theory  the  colour  of  a  patch  would  depend  on  the  colour  of 
the  neighbouring  patches,  this  is  consistent  with  Land’s  experiments.  Rubin  and  Richards 
(1984)  have  considered  ways  of  distinguishing  such  "material"  patches  by  considering  the 
spectral  properties  of  the  observed  intensities.  They  discuss  "lawful  processes"  (such  as 

1  The  rods  are  not  normally  assumed  to  contribute  to  colour  perception,  however,  we  will  arguo  later 
that  although  their  inputs  do  not  directly  contribute  to  perceived  colour  they  can  still  be  used  to 
"normalize"  colour  perception. 

’This  corresponds  to  the  standard  regularization  approach  (Poggio  and  Torre,  1984)  of  solving  an 
undercon3trained  problem  by  restricting  the  space  of  solutions. 
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orientation  and  shadowing)  which  change  the  albedo  while  still  corresponding  to  the  same 
material  and  propose  an  "opposite  sign  slope  condition"  to  detect  such  change. 

The  first  part  of  this  paper,  although  done  independently,  is  similar  to  work  of  Wandell  and 
Maloney  (1984a,  1984b)3.  In  their  work  the  illuminant  and  spectral  surface  reflectances 
are  also  expanded  in  terms  of  a  finite  number  of  basis  functions  in  spectral  space  and 
neighbouring  colour  patches,  with  the  same  assumed  illumination.  Our  formula  for  the 
number  of  patches  required  differs  from  theirs,  however,  and  furthermore  we  demonstrate 
a  simple  way  to  solve  the  equations.  Our  scheme  can  also  be  extended  using  results  from 
Hurlbert  and  Poggio(1985)  to  cases  when  the  spectral  components  of  the  illuminant  and  the 
reflectance  function  vary  with  position.  For  this  case  we  can  find  the  colour  of  an  object  at 
its  boundaries  and  then  fill  in  the  interior.  This  method  will  automatically  find  the  material 
boundaries  and  thus  will  act  as  a  material  edge  detector  in  the  sense  of  Rubin  and  Richards 
(1984). 

2.  The  Basis  Functions.  First  we  consider  the  case  when  there  is  no  spatial  dependence 
of  the  illuminant  or  the  surface  reflection.  The  inputs  to  the  receptors  are  then 

=  J  aM(X)E(X).9(X)r!X.  (2.1) 

We  now  expand  ft'(X)  and  9(X)  in  terms  of  basis  functions  7l,(X)  and  (7,(X) 


n  m 

ff(M  =  £  U  DM  S(X)  -  £  7>C,(X).  (2.2) 

«=1  J=l 

For  the  present  we  will  set  n  =  m  =  3  and  choose  the  same  basis  functions 

3  3 

£(X)  =  £  uDM  S(\)  ■=  £  7 jDM  (2.3) 

«=l  j—l 

This  is  merely  for  convenience  and  because  of  psychophysical  evidence  that  perceived 
colours  form  a  three  dimensional  space  (Helmholtz,  (1866),  for  a  review  see;  Feynman 
1965).  It  is  straightforward  to  generalize  the  analysis  to  different  values  of  n  and  m  and  for 
different  numbers  of  rods  and  cones,  Substituting  (2.3)  into  (2.1)  gives 

=  r,T,i  /«m(XW)«>(X)A  (2.4) 

where  we  use  the  summation  convention  over «  and  j. 

We  can  write  this  as 

3I  am  grateful  to  Ted  Adelson  for  bringing  their  work  to  my  attention. 


The  (T*)  depend  only  on  the  absorption  coefficients  and  the  basis  functions.  They  are 
therefore  parameters  of  the  system  and  hence  are  known  (the  system  might  be  taught  to 
learn  the  optimal  (?’£)).  The  (r()  and  (-yy)  describe  the  illuminant  and  the  surface  reflection 
respectively.  We  therefore  have  four  equations  (/<  ==  1  to  l)  for  six  unknowns  (»'  =  1  to  3 
,  j  =  1  to  3).  Note,  however,  that  there  are  effectively  only  five  unknowns  because  of  the 
scaling  ambiguity 

Ti  *-*  pr.-.Tfj  — •  (2.6) 

p 

This  ambiguity*  can  be  traced  back  to  (2.1)  where  it  corresponds  to  the  scaling 

E{\)  >-  PE{\),  S(\)  (2.7) 

P 

This  scaling  ambiguity  cannot  be  dealt  with  by  our  theory  and  a  further  normalization 
is  required.  Note  this  is  the  only  ambiguity  in  our  theory  and  only  one  normalization  is 
required.  This  is  in  contrast  to  the  theories  of  Land(1983),  Hom(1975)  and  Blake(1984) 
where  colour  is  computed  independently  in  separate  channels  each  of  which  has  to  be 
normalized  seperately.  Unlike  these  theories  (for  a  review  of  these  see  Hurlbert  1985)  our 
method  allows  interactions  between  channels  and  hence  we  only  need  one  normalization. 

To  resolve  this  ambiguity  we  can  normalize  the  light  source  by  making  (r.)  a  unit  vector,  or 
by  setting  n  =  1. 

A  single  patch  gives  us  four  equations  for  five  unknowns.  Suppose  we  have  a  neighbouring 
patch  with  a  different  spectral  surface  reflectance.  This  second  surface  reflectance  can  be 
expressed  in  terms  of  basis  functions  by 

3 

S*(X)  -  £  (2.8) 

<=l 

Now  if  we  make  the  reasonable  assumption  that  the  spectral  distribution  of  the  illuminant  is 
unchanged  on  the  second  patch  we  measure  four  new  quantities 

K  =  (2.9) 

4  whose  nature  was  clarified  during  discussions  with  Hurlbert  and  Poggio 
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which  involve  only  three  new  unknowns  (f>).  We  combine  this  with  (2.8)  to  get  eight 
equations  for  eight  unknowns.  Number  counting  suggests  this  is  enough  to  get  a  unique 
solution.  The  equations  are  non-linear  however  and,  as  we  will  show  in  the  next  section,  in 
some  cases  there  may  be  more  than  one  solution. 

For  the  human  eye  to  use  a  colour  scheme  like  this  with  four  receptors,  the  rods  are 
necessary  as  well  as  the  cones.  If  we  just  use  the  cones  and  have  three  receptors  it  follows 
from  (2.10)  that  we  can  only  have  two  basis  functions  and  hence  only  a  two  dimensional 
space  of  perceived  colour.  Observe  that  the  output  from  the  rods  can  be  treated  differently 
from  the  output  from  the  cones.  It  is  combined  with  them  to  solve  (2.4)  and  (2.6)  and 
determine  the  (f>)  and  the  (n)  but  it  is  only  the  inputs  from  the  cones  that  gives  the 
colours.  The  rods  and  the  cones  combine  to  determine  the  perception  of  colour  from  the 
spectral  information  at  the  cones. 

We  now  consider  the  general  case  when  there  are  p  receptors,  n  basis  functions  for  the 
illuminant  and  m  for  the  reflectance  function.  Suppose  we  have  q  patches.  Then  the  number 
of  equations  is  pq.  The  number  of  variables  is  n  +  mq  -  1  (the  l  comes  from  the  scaling 
ambiguity).  To  obtain  a  unique  solution3we  need 


pq  >  n  +  mq  —  I.  (2-10) 

This  is  a  neccessary  condition.  To  check  it  is  sufficient  we  will  show  that  the  non-linear 
equations  can  be  solved. 

Note  from  (2.10)  that  we  always  need  p  >  m,  that  is  the  number  of  different  photoreceptors 
must  be  greater  than  the  number  of  basis  functions  for  the  reflectance.  Thus  to  obtain  a 
three  dimensional  colour  space  with  this  scheme  at  least  four  photoreceptors  are  required.) 

The  number  of  patches  required  is 


q  > 


n  —  1 
p  —  m 


(2.11) 


So  if  there  are  two  basis  vectors  for  the  illuminant  and  surface  reflectances  (n  =  m  =  2)  and 
three  receptor  fields  (p  =  3)  we  only  need  one  patch. 


3.  Solving  the  Equations. 

We  now  consider  the  equations  in  more  detail.  They  are  of  the  form 


f>i  —  fiTfj  f'lj  I  k|i  —  TiijTij 


(3.1) 


5  This  equation  has  also  been  derived  by  J.  Rubin  (pers.  comm.).  II  does  not  imply  that  q  >  n,  a 
formula  stated  by  Wandell  and  Maloney  (1984a),  and  a  counterexample  is  given  in  the  next  section. 
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The  matrices  T{‘  are  parameters  of  the  system  given  by 

K  =  /  (3.2) 

If  the  basis  functions  (H,(\))  do  not  overlap  then  the  off-diagonal  terms  of  are  zero 
(Tf )  —  0,  i  ^  j).  This  case  is  considered  by  Huribert  and  Poggio  (1985).  This  simplifies  the 
equations  (3.1)  but  at  the  cost  of  preventing  interactions  between  different  colour  channels 
(although  interactions  can  be  introduced  later).  These  interactions  are  important  for  our 
purposes  since  they  mean  that  only  one  colour  normalization  (the  overall  scaling  given  by 
(2.7))  is  needed. 

The  basic  strategy  to  solve  the  equations  (3.1)  is  to  take  linear  combinations  of  the  matrices 
7’J‘  to  obtain  simpler  equations.  We  illustrate  this  by  considering  the  case  when  there  are 
three  photoreceptors,  and  two  basis  functions  («  —  m  -  2,  p  —  3).  As  shown  in  the  previous 
section  only  one  patch  is  needed  to  solve  the  equations.  These  are 

=  r,-7 W,  i  =  1,2.;  -  1, 2  M  =  l,  %  3.  (3.3) 

The  (T^)  are  three  two-by-two  symmetric  matrices.  The  space  of  such  matrices  is  three 
dimensional  and  so  the  form  a  basis  for  this  space.6  Hence  we  can  find  three  vectors 
p“  («  -  1,2,3)  such  that 


1  0 

.0  0. 


EM 


0  0 
.0  1. 


p!  T* 

,  *  t  J 


0  1 
l  0 


£««■ 


Combining  (3.4)  with  (3.3)  defines  three  (known)  numbers  A|,/t;,/1j  by 

iii 

=  E  —  tiTi»42  =  ^  Pm1*  =  fiTa  +r27l. 

M=t 

We  normalize  the  light  source  by  setting  rt  -- 1.  This  reduces  (3.5)  to 


(3.4) 


(3.5) 


_  7i  =  Ai.Tfjft  —  At,AiTt  -t-  7i  =  Aj.  (3.8) 

"We  choose  the  (T‘‘t )  so  that  they  are  linearly  independent.  We  are  free  to  do  this  since  they  are  the 
parameters  of  tho  system. 
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These  can  be  solved  to  give 


.  -4.1  ±  (/t3  - 

=  l,r2  - — ~  . - L — i h  =  -= 

2/1  / 


A3  *  (/t2  -  l/Ms)1/* 


(:t.7) 


This  gives  two  possible  solutions  (the  values  of  r2  and  ~i  -  are  interdependent).  The  constraint 
that  the  illuminant  and  the  reflectance  function  are  both  positive  functions  may  not  help  to 
resolve  this,  the  two  roots  of  both  r2  and  ~u  all  have  the  same  sign. 

Mote  that  only  one  patch  is  needed  to  solve  for  this  case  although  there  are  two  basis 
functions. 


The  General  Case 

The  basic  strategy  to  solve  the  equations  consists  of  taking  linear  combinations  of  the  ( T ?•) 
and  reducing  them  to  some  algebraic  equations  having  only  a  finite  number  of  solutions. 

For  four  photoreceptors  and  three  basis  functions  we  require  two  patches  with  the  same 
illumination.  We  have  two  sets  of  equations,  one  from  each  patch, 


=  rnjT^k.  =  r^T%.  (3.8) 

A  simple  way  to  solve  these  equations  is  by  choosing  the  basis  vectors  (8,{X))  so  that  the 
(7’5‘)  generate  a  space  spanned  by 


t 

0 

0 

0 

0 

O' 

0 

1 

O' 

0 

0 

0 

0 

0 

0 

t 

0 

1 

0 

> 

1 

0 

1 

) 

0 

0 

0 

.0 

0 

0. 

.0 

0 

0. 

.0 

1 

0. 

.0 

0 

1. 

The  conditions  on  the  basis  vectors  (/7,(X))  which  this  requires  are 

J  a*/J|(X)/Ja(X)fiX  =  0,  Ja'1//,(X)H2(X)«/X  =  J  a,lHt(\)lh{\)d\.  (3.10) 

To  see  this  observe  that  the  space  of  symmetric  three  by  three  matrices  is  six  dimensional. 
The  matrices  defined  by  (3.9)  define  a  four  dimensional  subspace  and  to  ensure  that  the 
(7’f;)  lie  in  this  subspace  we  must  impose  the  two  conditions  of  (3.10).  To  check  the 
exact  form  of  (3.10)  consider  all  possible  matrices  that  can  be  generated  by  linear 
combinations  of  (3.9).  First  observe  that  it  is  impossible  to  get  a  matrix  with  //J3  ^  0,  this 
gives  us  the  first  equation  of  (3.10).  The  second  equation  comes  from  noting  that  we  must 
have  II  |2  ~  llu  for  all  //. 

We  now  take  linear  combinations  of  (3.8)  using  four  vectors  </*  (a  —  1, 2, 3,1)  as  in  (3.4)  and 
(3.5)  until  we  generate  the  basis  matrices  of  (3.9).  The  equations  (3.8)  are  now  of  form 


Ai(r^r3)  +  A'i{n  +  r3)  +  /t3rij  -----  /\4t2t3)  /}|(r;r3)  +  B2(t3  +  t\)  +  /J3t|  =  n4r2r3, .  (3.13) 


The  (/!,)  and  (/>,)  are  linear  combinations  (using  the  (/“) of  the  (l,,)  and  the  (mM).  We  have 
normalized  by  setting  r,  =  1  and  divided  some  of  the  equations  by  r2  and  r3. 

This  gives 


_ (/!,  —  SDjTj 

(Ai  —  6Hi)t3  +  (/t3  —  SB3)  . 


(3.1-1) 


where  6  —  A2/H2 ■  We  can  substitute  this  into  (3.13)  to  obtain  a  quartic  equation  for  r3 
which  will  have  at  most  four  solutions.  Then  we  solve  (3.14)  for  r2  and  (3.11)  and  (3.12)  for 
the  (-7.)  and  (ft). 

The  basis  matrices  defined  by  (3.9)  were  choosen  arbitrarily.  The  same  analysis  could  be 
applied  for  other  basis  matrices.  For  some  choices  of  basis  there  is  only  one  possible 
solution.  An  example  is 
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This  basis,  however,  requires  that 


(MS) 


I  =  0  ,J  =  0.  (3.16) 

which  seems  an  unnatural  choice  for  the  (W,).  This  illustrates  the  important  point  that  the 
choice  of  basis  functions  is  crucial.  By  assuming  that  we  can  express  our  illumination 
and  surface  reflectance  in  terms  of  a  finite  set  of  basis  functions  we  are  performing  an 
implicit  type  of  regularization  (Poggio  and  Torre,  1984)  by  restricting  the  space  of  possible 
solutions.  Only  if  the  basis  functions  are  choosen  so  that  linear  combinations  of  them  are 
able  to  represent  most  of  the  naturally  occuring  illuminants  and  surface  reflectances  will 
methods  of  this  type  give  good  solutions. 
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This  strategy  of  taking  linear  combinations  of  the  (T?)  will  drastically  simplify  the  equations. 
It  is  also  possible  to  solve  equations  (3.8)  directly.  We  write  the  equation  as 


lM  —  +  T?2r2  +  T13t3)  +  'i-i(T3lrl  +  T22t2  +  T33t3)  +  '?3(T,3iri  +  T32t3  +  T33t3),h  —  1,  ..4 

(3.17) 

We  can  use  three  of  these  four  equations  to  solve  for  (-/*)  in  terms  of  the  (17).  This  will  leave 
us  with  two  polynomial  equations  relating  r2  and  t3  (we  have  set  u  =  1).  Similarly  we  solve 
similar  equations  for  (/?,)  in  terms  of  the  (a,)  leaving  another  polynomial  equation  in  r2  and 
r3.  These  two  equations  in  r2  and  r3  can  be  combined  to  give  solutions  for  r2  and  r3.  Since 
these  equations  are  polynomials  there  will  usually  be  more  than  one  solution,  and  in  some 
cases  an  infinite  number.  The  previous  method  involving  basis  matrices  is  a  lot  simpler.  If 
it  is  impossible  to  choose  a  basis  like  (3.9)  it  will  still  be  possible  to  simplify  the  equations 
by  taking  linear  combinations  of  the  (7'£)  to  reduce  the  matrices  to  simple  forms. 

Note  that  so  far  nothing  has  been  said  about  the  behaviour  of  the  solution.  It  is  even 
conceivable,  though  unlikely,  that  the  solutions  are  sometimes  not  positive  everywhere. 
They  are  guaranteed  to  be  smooth  provided  the  basis  functions  are.  Instead  of  choosing 
three  basis  functions  and  solving  for  their  coefficients  it  would  be  possible  to  have  more 
basis  functions  and  impose  some  a  priori  expectations  on  the  solutions.  Such  an  approach 
is  discussed  in  Richards,  Yuille  and  Poggio  (1985). 

The  Spatial  Extension 

In  this  section  we  consider  extensions  to  situations  when  the  surface  reflectance  and  the 
illumination  are  functions  of  space.  Hurlbert  and  Poggio  (1985)  have  shown  that  for  many 
cases  the  spatial  dependence  of  the  illumination  and  surface  reflectance  can  be  factored 
out  from  the  spectral  dependence  for  colour  patches.  For  these  cases  we  can  write 

S(x,X)  =  S(X)/(r)  (4.1) 

and 


/?(*,X)-=/-;(X)tfM  (4.2) 

where  J[x)  and  g(x)  depend  only  on  position. 

With  these  assumptions  the  spatial  dependence  of  the  illuminant  and  the  surface  reflectance, 
when  not  due  to  material  (colour)  changes,  is  factored  out  and  the  analysis  can  proceed 
as  in  the  previous  section.  There  is  an  important  difference,  however,  the  scaling  ambiguity 
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(2.7)  will  now  have  a  spatial  dependence.  In  other  words  we  will  not  be  able  to  distinguish 
the  spatial  dependences  /( z)  and  <j[x)  of  the  surface  reflectance  and  illumination.  For 
many  cases  the  illumination  will  have  little  spatial  variation  (for  example,  if  the  light  sources 
are  a  long  way  away)  and  we  can  set  g(x)  —  t  and  solve  for  f{x).  Otherwise  we  can 
consider  the  boundaries  between  two  materials  and  assume  that  the  spatial  dependence 
of  the  illumination  is  negligible  across  the  boundary.  We  can  then  solve  for  the  colour  on 
either  side  of  the  boundaries  (and  fill  in  if  necessary).  A  third  method  is  to  note  that  if  there 
is  a  single  light  source  the  spatial  dependence  <j(x)  will  have  a  definite  form  ,  typically  it  will 
fall  off  as  1  jr-  (whore  r  is  the  distance  to  the  light  source),  and  so  f(x)  can  be  determined. 

In  the  previous  sections  we  were  considering  a  Mondrian  style  world  consisting  of  flat 
patches  of  uniform  colour.  Our  method  would  sample  two  patches  and  obtain  two  equations 
(2.5)  and  (2.6).  For  the  more  general  case  we  would  sample  neighbouring  points.  If  the 
two  points  lie  on  the  same  material  the  l„  will  be  proportional  to  the  k „  and  the  equations 
will  be  underdetermined.  If  the  points  occur  at  a  material  boundary,  one  lying  in  each 
material,  then  we  have  enough  equations  to  solve  for  the  colours.  Since  we  are  considering 
neighbouring  points  it  is  reasonable  to  assume  that  the  illuminant  is  the  same  (including  the 
scaling  factor)  for  each  point. 

To  summarize:  we  sample  neighbouring  points  in  the  image  and  compare  their  values 
and  these  are  proportional  we  conclude  that  the  points  have  the  same  colour  (though 
not  necessarily  the  same  brightness)  and  hence  lie  on  the  same  material.  If  they  are  not 
then  we  are  on  a  material  boundary  and  can  solve  for  the  colours  on  either  sides  of  the 
boundary.  The  colour  cannot  change  inside  the  boundary,  or  else  we  would  detect  it,  and  so 
the  colour  on  the  boundaries  determines  the  colour  in  the  interior.  However  the  brightness, 
or  strength,  of  the  colour  is  undetermined.  If  the  illuminant  is  constant  in  space  we  can 
solve  this  by  a  single  global  normalization  (such  as  setting  rt  =  1,  as  in  (3.5)).  We  suggest 
making  this  assumption  even  if  the  illuminant  varies  spatially,  unless  the  variation  is  very 
large  or  there  is  strong  evidence  to  the  contrary  from  other  souces.  Thus  spatial  variation 
of  the  illuminant  will  be  interpreted  as  due  to  changes  in  the  reflectance  function  due  to 
orientation  and  other  changes. 

Other  ways  of  generalizing  to  allow  for  spatial  dependence  would  include  using  the  standard 
assumption  that  the  spatial  variation  of  the  illuminant  is  smaller  than  the  spatial  variation  of 
the  surface  reflectance  (Land  (1983),  Horn  (1974),  Blake  (1984)).  New  methods  for  doing 
this  are  discussed  in  Hurlbert  and  Poggio  (1985)  and  extensions  will  be  discussed  in  a 
forthcoming  paper. 
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